Collision of Viscoelastic Spheres: Compact Expressions 
for the Coefficient of Normal Restitution 
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The coefficient of restitution of colliding viscoelastic spheres is analytically known as a complete 
series expansion in terms of the impact velocity where all (infinitely many) coefficients are known. 
While beeing analytically exact, this result is not suitable for applications in efficient event-driven 
Molecular Dynamics (eMD) or Monte Carlo (MC) simulations. Based on the analytic result, here 
we derive expressions for the coefficient of restitution which allow for an application in efficient eMD 
and MC simulations of granular Systems. 

PACS numbers: 45.70.-n,45.50.Tn 



Introduction and description of the system. The collision 
of frictionless (smooth) viscoelastic spheres obeys New- 
ton's equation of motion, 



(«), 



(1) 



with the effective mass m e ff = m\mil (rri\ -\-m-i) and the 
compression £ = R\ + R2 — \rx — where r\ and r*2 are 
the time dependent positions of the spheres. F(. . . ) is 
the normal component of the vectorial interaction force 
F = F ■ e with the unit vector e = (r*i — F2) / \f\ — r*2|. 
For non-adhesive viscoelastic spheres it reads [l[ 

F = F cl + F dis = min (o, -pf /2 - ^Ap^i) (2) 



with 



3(l-i/ 2 ) 



(3) 



and Y, v and R c e stand for the Young modulus, the Pois- 
son ratio and the effective radius R c s = R1R2/ (R1+R2), 
respectively. The dissipative constant A is a function 
of the elastic and viscous material parameters [lj. The 
min(. . . ) function assures that the force is always repul- 
sive. 

The elastic part in Eq. (2), F , is the Hertz contact 
force Q while its dissipative part, _F dls , was first moti- 
vated in Q and then rigorously derived in [l|, Q , where 
only the approach in [l[ lead to an analytic expression of 
the material parameter A. 

While the knowledge of the interaction force, Eq. ([2]) 
is sufficient to perform Molecular Dynamics simulations 
(MD), the coefficient of restitution is needed to perform 
much more efficient event-driven MD and Direct Simu- 
lation Monte Carlo (DSMC) as well as for the Kinetic 
Theory. By disregarding the dynamics of the collision 
process and idealizing the collision as an instantaneous 
event, the coefficient of restitution relates the postcolli- 
sional normal velocity, to the normal component of 



the (precollisional) impact velocity, v, 
e=-i'/v. 



(4) 



In general, the coefficient of restitution is not a con- 
stant but depends on the details of the interaction force 
and the impact velocity. It can be obtained by integrat- 
ing Eq. ((T|) with the initial conditions £(0) = and 
£(0) = v, assuming that the spheres start contacting at 
t = 0. The coefficient of restitution is then obtained from 



<(tc)/v : 



(5) 



where the duration of the collision, t c , is determined by 
the condition 



|(i c ) - t c > , 



(0) 



that is, the collision terminates at time t c when the in- 
teraction force vanishes. 

Solving the set of equations (|5|6|) is a complicated prob- 
lem which was solved rigorously in [5| . The solution reads 



e = 1 + J> (> / V/ 10 ) fe ^l + ]T h k vl (7) 
where we define the shorthand v* and with 
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(8) 



This solution is exact since all coefficients hk are analyt- 
ically known (see [Ij]). It is, moreover, universal since all 
material and particle properties are covered by j3, that 
is, the hk are pure numbers which are independent of the 
material and particle properties. 

Albeit exact, there are two main problems with the so- 
lution, Eq. ([7]), which prohibit its application in efficient 
MD or DSMC simulations: First, it converges extremely 
slowly. To obtain e up to quadratic order in v we need 
20 terms of the series expansion. Second, wherever we 



2 



0.8 



0.6 



0.4 



0.2 



// .' 
' i 

i 



k=7 

k c =8 

k=9 

k c =10 

k c =I5 

k=20 

[1/4] e 

numeric (exact) 



✓ / 



\ s 



\ ^ \ N 



0.25 



0.5 0.75 

v. 



1 



1.25 



FIG. 1. (color online) Coefficient of restitution, e, over 
v» = [3 1/2 v 1/w . The analytic solution, Eq. 0, truncated at 
different order k c leads to divergence. The dotted line shows 
e as it follows from the numerical solution of Eqs. (|5I6|I . It 
almost coincides with the thick green line showing the Pade 
approximant [l/4] e , Eq. to the analytical solution, Eq. 

(for explanation see the text below). 



solution is exact and may serve as a benchmark for our 
approximative solution, even for large values of . Using 
the numerical solution we find the asymptotical behavior 



lim e(v*) = v 

v v — s-oo 



-3.2 



(10) 
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truncate the series at some order k c , Eq. (J7]) diverges to 
e — > ±oo, depending on the sign of h^. 

The divergence of the truncated series is a serious 
problem: Given the very accurate experimental data by 
Bridges et al. 6] for the coefficient of restitution of ice 
balls at very low temperature whose material and par- 
ticle properties correspond to 1.307 (sec/cm) 1 ^ 5 . From 
Fig. [T]we see that the series truncated at order k c — 20 
starts deviating at v„ w 1 corresponding to the impact 
velocity v = v l°//3 5 « 0.262 cm/sec. That is, for typical 
impact velocities of v ~ 1 m/sec we would need to go to 
impractical high truncation order. 

From an approximative expression for the coefficient 
of restitution for applications in efficient MD and DSMC 
simulations, we request that a) the approximative solu- 
tion is close to the correct solution, b) it can be computed 
efficiently, that is, it contains only a small number of uni- 
versal coefficients which are independent of the material 
and particle properties, and c) the representation must 
not reveal divergencies unlike the truncated series, Eq. 
0, shown in Fig. [TJ 

Numerical solution. As described in 0,0, Eq. Q with 
the interaction force Eq. ([2]) and the corresponding initial 
conditions may be scaled to 

x + x 3/2 +vl±y/x = 0, x(Q) = 0, ar(0) = l (9) 

with the only free parameter v* = /3 1 / 2 ?; 1 / 10 . Compres- 
sion and time are scaled by x = ^/[(p/m cS )- 2 ^v^ 5 } 
and t = t/[(p/m e fi)~ 2 / 5 v~ 1 / 5 ]. From the numerical so- 
lution of Eq. ([9]) we determine e(v*) via Eq. ([5]): e — 
—£(t c )/v — —x{t c ), where r c is obtained from the condi- 
tion x(t c ) — 0, t c > 0. Apart from numerical errors, this 



FIG. 2. (color online) Coefficient of restitution e for large 
v* . The thick green line shows the numerical solution of Eq. 
revealing the asymptotic behavior e = u^ 3 ' 2 (dotted line). 
Additionally various Pade approximants, Eq. (| 1 1[) . of the 
analytical solution, Eq. are shown (discussion see text 
below). The Pade approximants [3/6] e and [15/18] s (virtually 
identical) agree almost perfectly with the exact solution. 

Pade approximants. Using the analytical solution, Eq. 
((7|), and the asymptotics, Eq. ([TO)) , we construct an ap- 
proximative expression for s(v) which agrees with the 
analytical solution for the entire range of definition, 
v € (0, oo), and is thus much more suitable for numerical 
simulations. The Pade approximant [m/n] e (v*) approxi- 
mates the m + n times differentiable function e(u*) by a 
rational function 

[m/n] e („.) = ^r°H (11) 

in a way that the Maclaurin series of the approximant 
and of the approximated function match up to order 
to + n: e(0) = [m/n] £ (0), e'(0) = [m/n]' e (0), 
e (™+«)(0) = [TO/n]i m+n) (0). Asymptotically, the Pade 
approximant behaves like a power law, lim 1Ii> _ i . 00 [TO/?i] e ~ 
v™~ n . These properties allow to represent the function 
e(««) similar to a Taylor expansion for small arguments 
and asymptotically as a power law, thus, convergent if 
to < n, see Ref. @. 

Since e ~ v° with a m -3 (see Eq. JTO]) and Fig. [2]) 
we chose a Pade approximation [to/to + 3] e . To find an 
accurate yet compact approximant to Eq. ([TJ we start at 
to = and increase the order until sufficient agreement 
with the exact solution is achieved. The result is shown 
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in Fig. [2j [0/3] e is certainly not acceptable, [l/4] e offers 
a good tradeoff between simplicity and accuracy. [2/5] s 
reveals a pole at v* « 5.68, therefore, it is suitable only 
for small impact velocity, < 10 3 . For ice spheres as 
described in Q this implies v < 2.6m/sec. The next 
order, [3/6] e , offers almost perfect agreement with the 
benchmark. We checked all orders up to [25/28] e and 
could not find any significant improvement as compared 
to [3/6] e . As an example, [15/18] e is shown in Fig. [2] 

Table Q] displays the coefficients aj and bi for the rel- 
evant Pade approximants [m/m + 3] e , m G {0,1,2,3} 
and Fig. |3] shows these Pade approximants together with 
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TABLE I. Coefficients of the Pade approximants [m/m + 3] e 
for m £ {0, 1, 2, 3}. [2/5] e reveals a pole at v* ~ 5.6801. 

the exact (numerical) solution. Again [1/4] e and [3/6] e 
turn out to be good compromizes between accuracy and 




FIG. 3. (color online) Coefficient of restitution e over v*. 
The first four Pade approximants are shown together with the 
numerical (exact) solution. The inset shows a magnification. 
The order [3/6] e (dotted line) coincides almost perfectly with 
the exact solution in the entire range of definition. 



Conclusion. The universal exact solution, Eq. (J7J, 
for the coefficient of restitution of smooth viscoelastic 
spheres cannot be applied directly in eMD and DSMC 
simulations since the series diverges for any finite trun- 
cation order. We have shown that the Pade approxima- 
tions of order [1/4] e and [3/6] e are suitable to represent 
the coefficient of restitution over the entire range of im- 
pact velocities including its asymptotic behavior up to 
an excellent accuracy and we provided the constants of 
this approximation. Similar as the full solution, Eq. (|7|). 
the Pade expansion is universal, that is, the constants 
di and bi are universal. They neither depend on mate- 
rial properties (Young modulus, Poisson ratio, dissipa- 
tive constant) nor on particle properties (radii, masses). 
All non-universal parameters enter exclusively via (3, Eq. 
([8j), which in turn enters the argument of the Pade expan- 
sion via v* = /3 1 / 2 !; 1 / 10 with v being the impact velocity 
in physical units (cm/sec). Thus, the presented Pade 
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FIG. 4. (color online) Coefficient of restitution e as a function 
of the impact velocity v. The Pade approximant [3/6] e (dot- 
ted line) agrees almost perfectly with the numerical integra- 
tion of Newton's equation, Eqs. (JTJ6]) in the entire range of im- 
pact velocity, v (physical units), while the analytical solution, 
Eq. (0, truncated at order as large as k c = 40 diverges at 
v w 0.3 cm/sec. For the material constant, 1.307 (sec/cm) 1 '' 5 , 
we used the experimental values by Bridges et al. [y] for the 
collision of ice spheres at low temperature. 

The precision of the approximant can be assessed in 
Fig. @] which shows the Pade approximation together 
with the numerical integration of Newton's equation, Eq. 
(TJJ) in combination with Eqs. dUS]), and with the diver- 
gent analytical solution, Eq. (0, truncated at order as 
large as k c = 40. We see that over the entire range of 
definition, the Pade approximation coincides almost per- 
fectly with the numerical solution and with the truncated 
analytical solution up to w « 0.3 cm/sec where it starts 
to diverge. For the material constant, 1.307 (sec/cm) 1//5 , 
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we used the experimental values by Bridges et al. [6| for 
the collision of ice spheres at low temperature. The corre- 
sponding data is also shown in the plot. While the agree- 
ment between the exact analytical result, the numerical 
integration and the Pade approximant is remarkable, the 
experimental data slightly deviates. This deviation is 
not surprising since besides viscoelasticity, described by 
the force Eq. ((2|), other forces may contribute, such as 
surface forces, plastic deformation, adhesion etc. 
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